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Background: This paper describes a 2D and 3D simula- 
tion engine that quantitatively models the statics, dynamics, 
and non-linear deformation of heterogeneous soft bodies in 
a computationally efficient manner. There is a large body of 
work simulating compliant mechanisms. These normally as- 
sume small deformations with homogeneous material prop- 
erties actuated with external forces. There is also a large 
body of research on physically -based deformable objects for 
applications in computer graphics with the purpose of gener- 
ating realistic appearances at the expense of accuracy. 
Method of Approach: Here we present a simulation frame- 
work in which an object may be composed of any number of 
interspersed materials with varying properties (stiffness, den- 
sity, etc.) to enable true heterogeneous multi-material simu- 
lation. Collisions are handled to prevent self-penetration due 
to large deformation, which also allows multiple bodies to 
interact. A volumetric actuation method is implemented to 
impart motion to the structures which opens the door to the 
design of novel structures and mechanisms. 
Results: The simulator was implemented efficiently such 
that objects with thousands of degrees of freedom can be 
simulated at suitable framerates for user interaction using a 
single thread of a typical desktop computer. The code is writ- 
ten in platform agnostic C++ and is fully open source. 
Conclusions: This research opens the door to the dynamic 
simulation of freeform 3D multi-material mechanisms and 
objects in a manner suitable for design automation. 
Keywords: Soft body, non-linear mechanism, heterogeneous 
materials 



1 Introduction 

Scientific physics simulators are traditionally used to 
model small deformations of homogeneous linear-elastic 



materials. Recently, multi-material additive manufacturing 
methods have been developed that fabricate heterogeneous 
objects out of two or more materials. The properties of 
these co-fabricated materials can range from rigid plastics 
and metals to very soft rubber with linear deformation greater 
than 200% (TTJ. The inclusion of these soft, rubbery mate- 
rials necessitates the consideration of large, non-linear ge- 
ometric deformations to accurately predict the physical be- 
havior. Because hard and soft materials can be internally 
combined and patterned in 3D with very few constraints, a 
new paradigm of physics simulation becomes necessary to 
efficiently predict the combined material properties and dy- 
namic behavior. 

There are many established methods and implementa- 
tions for simulating deformable soft bodies (8). Consider- 
ing the large deformations and relatively low stiffness of the 
materials involved, the physically-based dynamics are often 
significant and must be modeled (Figure [T}. Much of the 
development in simulating soft bodies has been driven by 
the computer graphics community. Many of the well es- 
tablished physics engines provide support for dynamic de- 
formable bodies, whether ID rope, 2D cloth, or 3D "jello" 

cm). 

The goal of these simulations is generally to create re- 
alistic visual effects in real time at the expense of accu- 
racy (18]|2||7]|3). For instance, lattice shape matching (15) 
creates visually appealing rubber effects very efficiently and 
is unconditionally stable. However, the underlying methods 
are geometrically based, which limits their direct application 
to quantitative engineering analysis problems. Other simu- 
lators are derived from more physically based principles |5 ], 
but their performance at predicting real-world behaviors is 
unverified. Deformable body simulators have also been de- 
veloped specifically for real-time surgery simulation 1 1 , 20|. 




(a) (b) 

Fig. 1 . Two frames of a hollow ball bouncing under gravity illustrate 
the undeformed ball just before impact (a) and the highly deformed 
ball the moment of highest deformation (b). 



These simulators address challenges such as modifying the 
geometry dynamically to simulate incisions, but due to the 
variance of biological materials it is also difficult to verify 
quantitative accuracy. 

There is also a large body of work regarding the sim- 
ulation of compliant mechanisms and design thereof (6|[T2). 
Existing efforts focus on small displacements 1 10 ] or discrete 
thin beam members that can flex significantly |16||4). In the 
simulation framework presented here, an entire freeform 3D 
shape can be non-linearly deformed, leading to many novel 
possibilities for soft mechanisms that cannot be simulated 
efficiently using current techniques. 

1.1 Finite Element vs. Mass-Spring methods 

Finite element analysis (FEA) is a well established 
method of simulating the mechanical behavior of objects. 
Advantages include the ability to solve a system with irreg- 
ularly spaced discretized mesh elements. A stiffness matrix 
is composed containing information about the connectivity 
of the entire mesh and the local material properties at each 
node. However, this system can only be efficiently solved 
if the underlying equations are linear. Thus, deformations 
that change the geometry significantly require periodic re- 
meshing [20]. Other non-linearities such as friction and ad- 
vanced material models require additional levels of iteration 
to solve. 

Mass-spring methods are widely used for deformable 
bodies, especially in dynamic simulations for computer 
graphics (8) . Advantages include relative simplicity and han- 
dling large deformations and other non-linearities with ease. 
An object is decomposed into discrete point masses con- 
nected by springs. Thus, the entire system forms a system 
of ordinary differential equations (ODEs) that can be inte- 
grated directly to solve for the behavior of the system. This 
makes these particle-based physical simulations very com- 
putationally efficient at the expense of accuracy. 

1.2 Freeform Mesh vs. Voxels on a Lattice 

There are a number of tradeoffs associated with choos- 
ing either a freeform mesh or a lattice of voxels to dynam- 
ically simulate a heterogeneous object. Both FEA methods 



and many existing gaming physics simulators use a freeform 
mesh to discretize a 3D object for simulation. By allowing 
the vertices to lie at any position within the object, there is 
greater control over the local detail of the simulation. Specif- 
ically, this allows objects to be meshed based on the desired 
accuracy in a given region or dynamically re-meshed based 
on the current regions of interest in a deformed shape (9). 
Care must be taken when forming the mesh such that the as- 
pect ratio of each element does not vary significantly in order 
to preserve accuracy. However the advantages of freeform 
meshes quickly diminish as materials of different stiffnesses 
and properties are interspersed within the object. This con- 
strains the mesh generation process and can potentially cre- 
ate very large and inefficient meshes, such as the case of a 
dithering between two materials. 

Limiting the discretized elements in a simulation to vox- 
els has a number of favorable advantages. This approach en- 
ables efficient computation of the force of each constituent 
element, since they begin on a principle axis with identical 
lengths. Additionally, the stiffness of each linking beam can 
be pre-computed based on the stiffnesses of each constituent 
voxel so that each individual voxel can have a unique stiff- 
ness without altering the efficiency of the simulator (Figure 
[2|. This allows heterogeneous materials to be simulated with 
the same computational complexity as homogeneous materi- 
als. Additionally, using a voxel lattice eliminates the possi- 
bility of ill-formed meshes 1 17 ]. However voxel lattices are 
at a disadvantage to freeform meshes when large regions of 
homogeneous material are present or very fine local details 
must be simulated. 

1.3 Application in Design Automation 

Because of the exponentially increasing design space 
enabled by multi-material additive fabrication methods, de- 
sign automation will play an increasing role in the design and 
optimization of structures that fully take advantage of these 
capabilities. Most design automation algorithms depend on 
many, many physical evaluations. Therefore a balance must 
be struck between calculating accurate results while mini- 
mizing CPU cycles. In the simulation framework presented 
here, static and dynamic properties are quantitatively very 
close to the analytical solutions for simple textbook scenar- 
ios. Features such as collision detection are in place to avoid 
the great inaccuracy of self intersection, but are not meant 
to draw scientific conclusions about the interaction between 
two soft bodies. By carefully budgeting CPU cycles, the 
simulator can accurately model physical properties while not 
wasting undue time on negligible effects. 



2 Physics Engine 

2.1 Heterogeneous Deformable Body Core Simulator 

A voxel-based mass-spring lattice was chosen to best 
simulate the dynamics of highly deformable heterogeneous 
materials. Several measures were incorporated to mitigate 
the lack of quantitative accuracy normally associated with 
the discrete particle-based simulations. Each lattice point 
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Fig. 3. Each voxel is modeled as a lattice point with mass and ro- 
tational inertia (red). Voxels are connected by beam elements (blue) 
with appropriate translational and rotational stiffnesses leading to re- 
alistic deformation under applied forces and moments. 

implies that even though the physics engine presented here 
is capable of modeling large aggregate non-linear deforma- 
tions, the accuracy drops off as the angle between any two 
adjacent voxels becomes too large for a reasonable small an- 
gle approximation. 



(b) 



(c) 



Fig. 2. Advantages of a lattice-based voxel simulation include the 
native ability to simulate objects with multiple interspersed materials 
of varying properties (a). Here, the blue material is 100 times stiffer 
than the red, leading to higher deformation along the top (b). Internal 
deformation (strain) is also shown (c). 



was modeled with six degrees of freedom. In addition to the 
traditional three translational degrees of freedom, all three 
rotational degrees of freedom are stored and updates as part 
of the state. 

It follows that not only is a mass stored for each voxel, 
but also its rotational equivalent moment of inertia. Instead 
of using simple extension springs to connect adjacent points, 
more complex beam elements were used that resist lateral 
shearing and rotation in all axes in addition to extension. By 
setting the properties of the beam equal to the equivalent size 
and stiffness of the bulk material connecting two voxels, a 
good approximation of the aggregate bulk material behavior 
is obtained. 

To prepare a given geometry for simulation, the target 
geometry was first voxelized into cubic voxels. Each subse- 
quent relaxation step consists of two steps: (1) calculating 
all internal forces, then (2) updating all positions. In order to 
preserve the proper dynamics of the system positions were 
updated synchronously so that the order of calculation is ir- 
relevant. In order to capture information about both transla- 
tion and rotation of the voxels relative to each other, a con- 
stant cross-section beam element was used to connect adja- 
cent voxels in the lattice. Beam elements resist both trans- 
lation and rotation by exerting biaxial bending, transverse 
shear, and axial stretching forces in response to appropriate 
displacements (Figure [3]). Here we use a standard Bernoulli- 
Euler beam theory. It is important to note that the Bernoulli- 
Euler beam theory assumes a linearized beam model. This 



2.1.1 Compositing Adjacent Dissimilar Materials 

The Bernoulli-Euler beam theory also requires the mate- 
rial to be elastic and isotropic. When two adjacent voxels are 
composed of the same material, the elastic modulus and stiff- 
ness of this material are used in the equivalent beam connec- 
tion. However, when the materials have differing properties, 
an appropriate composite property must be calculated. To 
this end we approximate the composite stiffness of a bond 
between two dissimilar materials by: 



E r = 



2E X E 2 
E { +E 2 



(1) 



where E c is the composite elastic modulus and E\ and E 2 
are the two constituent stiffnesses. Since the elastic modulus 
directly corresponds to the spring constant of each bond, this 
is analogous to combining two springs of half the length in 
series with dissimilar stiffnesses. The composite shear mod- 
ulus are calculated in a similar manner: 



2GiG 2 
Gi+G 2 



(2) 



where G c is the composite shear modulus and G\ and 
G2 are the two constituent shear moduli. Because Poisson's 
ratio relates the elastic modulus to shear modulus, it follows 
that the composite Poisson's ratio ju c is calculated by: 



E^ 

2G r 



(3) 



2.1.2 Simulation elements 

Since solid objects are represented by a network of 
beams connecting nodes, the physical parameters for these 



beams must be calculated. Because the geometry is con- 
strained to voxels, the length / of the beam was taken to be 
the distance between the voxels, and the cross-sectional area 
of the beam A was I 2 . The standard formula was used to cal- 
culate the bending moment of inertia (7), given that in this 
case both the base b and height h equal the lattice dimension 
/. 



bh 3 
~V2 



12 



(4) 



The torsion constant (7) was approximated by the polar 
moment of inertia of a rectangular cross-section beam, and 
calculated as 



bh(bb + hh) 
12 



6 



(5) 



Using the standard Hermitian cubic shape functions for 
beam elements the stiffness matrix was determined for a 
beam element with 12 degrees of freedom: Three transla- 
tional and three rotational degrees of freedom for each end- 
point of the beam. This can be assembled into a stiffness 
matrix. The result for a beam element oriented in the posi- 
tive X direction is as follows: 
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In order to reduce the computational complexity, the 
transformations to rotate each individual element from its 
initial orientation into the positive X direction were precom- 
puted. Due to the cubic voxel lattice constraints, all elements 
are initially located on principal axes. This makes the trans- 
formation into the X-direction stiffness matrix computation- 
ally trivial. As the structure deforms over the course of a 
simulation, an additional transformation was calculated and 
applied to each element to translate the position and angle of 
the first voxel to zero. Therefore X\,Y\,Z\, Q Xl , Q yi , and Z] 
all become zero and drop out of the calculation. The large 
matrix calculation above is reduced to: 



F Xl = -a\X 2 

F yi =-hY 2 +b 2 a Z2 

F ZI = -biZ 2 +b 2 Q y2 
Mq xi = -a 2 Q X2 
Mo n =b 2 Z 2 +b 3 Q n 
Mq zi = -b 2 Y 2 + b 3 Q Z2 
F X2 = ~E Xl 
F y2 = —F yi 



m$ X2 = a 2 Q X2 
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M Qz2 =-b 2 Y 2 + 2b 3 e z2 
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where the stiffness matrix [K] is 
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The resulting forces and moments are then transformed 
back to the current orientation of the bond using the inverse 
of the transform calculated to arrange them in the positive X 
axis. The forces for all the bonds are calculated separately, 
then total forces (F t ) and moments (M t ) on each voxel are 
summed according to how many bonds n are connected to it. 



M t 



where 



b=n 

b=0 
b=n 

b=0 



(25) 



(26) 



2.1.3 Integration 

Because momentum plays a key role in all dynamic sim- 
ulations, two integrations are necessary to update the posi- 
tion realistically. For this physics engine, double Euler inte- 
gration was used. Although there are more accurate integra- 
tion methods, such as the Runge-Kutta (RK4) method, Euler 
was chosen because the massive number of discrete voxels 
and non-linear effects such as stick- slip friction are not well 
suited for the predictive steps of the RK4 integration scheme. 
The state of each voxel was represented by three dimensional 
position (D) and rotation (0) vectors and three dimensional 
linear (P) and angular momentum ($) vectors. In order to 
advance the simulation from time t n to time t n +\ =t + dt, 

Pt n+l =Pt„+F t dt (27) 
D tn+1 =D tn + ^dt (28) 
l n+l =ty tn +M t dt (29) 

where m is the mass of the voxel and / is the rotational 
inertia. 

2.1.4 Choosing the Timestep 

A critical aspect of implementing a robust physics sim- 
ulation driven by Euler integration is to choose a suitably 
small timestep to prevent numerical instability. However, in 
order to be computationally efficient, the timestep should not 
be unduly small. Fortunately, it is trivial both conceptually 
and computationally to determine the longest stable timestep 
at each iteration of the simulation. In an oscillating system, 
the simulation will be stable if 



dt<—*— (31) 

Because each bond between voxels is essentially a mass- 
spring-damper system, C0o m is simply the maximum natural 
frequency of any bond in the system. The stiffness of each 
bond was divided by the minimum mass of either voxel con- 
nected to it to calculate the maximum natural frequency of 
each bond according to 




where kt, is the stiffness of the bond and m m is the mini- 
mum of either mass connected by this bond. 



2.1.5 Damping 

Once an optimal timestep has been chosen, it is neces- 
sary to implement damping into the system to avoid the accu- 
mulation of numerical error as well as to enable realistically 
damped material properties. Because one application goal of 
this simulation involves unconstrained motion of soft bodies, 
damping must be included at the local interaction between 
voxels, not just applying a force proportional to each voxel's 
global velocity which would damp rigid body motion. 

The local damping between adjacent voxels ensures that 
modal resonances at the scale of a single voxel do not ac- 
cumulate. For each bond between two voxels, A force was 
applied to each voxel opposing the relative velocity between 
them. However, because rotational degrees of freedom al- 
low this bond to be spinning, both angular and translational 
velocities must be correctly accomodated to make sure rigid 
body motion is not being damped. For each bond, first the 
average position, velocity, and angular velocity were calcu- 
lated. Then the velocity of the second voxel relative to the 
first (V2->i) is calculated according to 

V 2 ^l = (Vl-Va) + 02 -D a ) X CO. (33) 

where V2 is the velocity of the second voxel, V a is the 
average velocity of two voxels, D2 is the position of the sec- 
ond voxel, D a is the average position, and CO. is the average 
angular velocity. This is in effect subtracting out the rigid 
motion components of the relative velocity such that they are 
not damped. Then for each each voxel the damping force is 
calculated according to the standard linear damping formula: 

F d = l^VmkVy (34) 

where F d is the damping force to be applied to the voxels 
with mass m attached to a bond with stiffness k and a rela- 
tive velocity V r . The damping ratio £ is normally selected to 
be 1, corresponding to critical damping. Likewise, angular 
velocities are also damped according to 

M d = 2^^Ik^ r (35) 

where M d is the damping force to be applied to the vox- 
els with a rotational moment of inertia / attached to a bond 
with rotational stiffness k and a relative angular velocity G) r . 
The rotational damping ratio £ is also normally selected to be 
unity. However, even though each bond is critically damped 
locally, the structure as a whole is still quite underdamped. 
So, each voxel was also variably damped in a similar manner 
relative to ground according to the situation at hand. 

2.2 Collisions 

2.2.1 Gravity, floor and friction model 

In order to properly simulate freely moving soft bod- 
ies, gravity is necessary. As the force is summed on each 



voxel, the mass of the voxel times the acceleration of gravity 
was subtracted from the vertical component of force. In con- 
junction with gravity, a floor was implemented to for objects 
to rest on. Because the maximum simulation time step that 
can be taken is limited by the maximum stiffness between 
any two connected masses, the effective normal stiffness of 
the floor on any voxels in contact with it cannot be infinitely 
high. In order to keep the simulation as efficient as possible, 
the stiffness of each voxel contacting the floor was the stiff- 
ness of the floor in that location. Although this allows signif- 
icant floor penetration in some cases, the qualitative behavior 
is appropriate. Potential collisions with the floor are trivial 
to detect by simply comparing the vertical position of each 
voxel to the ground plane, after accounting for the current 
size of the voxel. 

Although a standard linear friction model would provide 
a relatively realistic simulation, much more interesting and 
realistic behavior can be observed using a Coulomb friction 
model. This implies that a voxel at rest with the floor will 
resist any motion until 



\Fi\>ju s F n (36) 

where F\ is the horizontal force parallel to the ground, ju s 
is the coefficient of static friction between the voxel and the 
ground plane and F n is the normal force pressing this voxel 
into the plane of the ground. A boolean flag is set indicating 
to the simulation that this voxel should not move laterally, 
but can still move in the direction normal to ground such 
that it can be unweighted and then moved laterally. Once the 
static friction threshold has been exceeded at any given time 
step, the voxel is allowed to begin motion in the appropriate 
lateral direction by clearing the boolean static friction flag. 
The voxel is allowed to move in all three dimensions, but 
a friction force is applied opposing the lateral direction of 
motion according to: 



\Fi\=ii d F n (37) 

where ju^ is the dynamic coefficient of friction. In order 
to properly detect when a voxel has stopped lateral motion, a 
minimum motion threshold must be set. Otherwise the voxel 
will never re-enter the static friction state until the velocity is 
less than the precision of a floating point variable. In order 
to detect a stopping voxel, especially one that would change 
direction and incorrectly bypass the effects of static friction, 
a voxel is artificially halted if 



Vl < (38) 
m 

where m is the mass of the voxel in question. Because 
the force of friction (F n jUd) is always directly opposed to the 



voxels lateral velocity (V/), the voxel is stopped if the pro- 
jected change in velocity would change the direction of the 
surface velocity, which would involve the momentary stop- 
ping of the voxel. Collisions are also damped normal to the 
direction of contact with a user variable damping ratio rang- 
ing from zero (no damping) to 1 (critical damping). 

2.2.2 Self collision detection and handling 

Collision detection between voxels must be imple- 
mented carefully to avoid this being a bottleneck in CPU 
cycles. Especially in simulations with many independently 
moving particles, the 0(n 2 ) process of checking every par- 
ticle against every other to detect collisions is prohibitively 
expensive in CPU cycles. In a voxel simulation such as this 
there are many ways we can make collision detection more 
efficient. Since large deformations and multiple bodies are 
possible, we cannot simply exclude collision detecting be- 
tween voxels that are connected. However, we can imme- 
diately assume that voxels on the interior of an object my 
be disregarded for any collisions, assuming (as we will) that 
collisions are handled in such a way that overlaps cannot pen- 
etrate the outer shell. Upon import into the simulation, a list 
of surface voxels is precomputed, since this information will 
never change. 

The next step is to build a list of voxel pairs that are 
within a collision horizon (reasonable range) of each other. 
Voxels on this shortlist should be compared at every timestep 
for potential overlap. The collision horizon was chosen to be 
a distance equivalent to two voxels. However, it is unde- 
sirable to watch for potential collisions between voxels that 
are adjacent and connected in the lattice, since the internal 
forces between them already resist penetration. To account 
for this, a list of voxels within a 3D manhattan distance of 3 
in the lattice is precompiled upon import into the simulation 
for each voxel. Once complete, any potential collision in- 
teractions can be compared against this list in linear time to 
exclude the computational overhead of calculating spurious 
collision interactions. 

Finally, the list of potential voxel collision pairs must 
be updated often enough that out-of-range voxels beyond the 
collision horizon do not have a chance to penetrate before be- 
ing recognized as potential collisions. To accomplish this, a 
global maximum motion variable is initialized in the simula- 
tion. Each time step, the magnitude of the maximum velocity 
of any voxel in the simulation is added to this variable. While 
the maximum motion is less than half the collision horizon, it 
is guaranteed that no voxels can overlap into a collision. This 
is extremely conservative, but also computationally trivial to 
compute. 

2.3 Volumetric Actuation 

For convenience, we will refer to volumetric actuation in 
the context of materials with a non-zero coefficient of ther- 
mal expansion (CTE) in conjunction with a changing "tem- 
perature" control variable. However, volumetric actuation 
may be physically achieved in a variety of ways, so there is 
no reason to assume that the results presented here are appli- 



cable only to temperature changes. There is also no reason 
that different materials within the simulation couldn't expand 
or contract out of sync to multiple independent control vari- 
ables. This would be analogous to having multiple "tempera- 
tures" that only affect certain materials. But in the following 
discussion, we will refer to only a single temperature control 
variable. 

With the soft body relaxation engine in place, such vol- 
umetric actuation is implemented by simple changing the 
nominal rest length between adjacent voxels when comput- 
ing the elastic force between them. If the elastic force (Fe) 
between two voxels is normally calculated according to 



F e =k(P 2 -P 1 -D Npi ^ (39) 
to add in the effects of volumetric actuation, 



D NPl ^ |r=r c =(l + ^4^(r c -r r ))D^^ 2 \ T=Tr (40) 



where L>np 1 ^p 2 \t=t c is the modified rest distance based 
on the current temperature, OCi and OC2 are the coefficients of 
thermal expansion of the bond's constituent materials, T c is 
the current temperature and T r is the reference temperature, 
at which there is no temperature-based expansion or contrac- 
tion. 



3 Validation 

In order to ensure that the physics engine was perform- 
ing properly, we compared both static and dynamic behaviors 
of cantilever beams to finite element and analytical solutions. 
To verify the static behavior of the simulation, beam deflec- 
tions of both thin and think cantilever beams were compared 
to a linear direct stiffness method and (in the case of the thin 
beam) to the analytical solution. The results are outlined in 
Table [T] For the thin beams, 20x1x1 voxels were used with 
a physical size of 1mm each, for a total beam size of 20mm 
long by 1mm thick. A material stiffness of IMPa was spec- 
ified. The force at the end of the beam was selected to be 
0.03mN so that the displacement would be small (less than a 
voxel-height). This ensures that small angle approximations 
of the analytical solution are valid. The thick beams were 
modeled as 10x5x5 blocks of 1mm voxels with the same 
stiffness. In this case 0.1N of force was applied at the free 
end to achieve a non-negligible displacement. 

The non-linear mass/spring method presented here re- 
sults in slightly smaller displacements than the direct stiff- 
ness method and the analytical solutions (Table [T}. This dif- 
ference is negligible in the thin beam case. The difference 
is more pronounced in the thick beam case. This is likely 
because the deformation (Figure [4} is large enough that the 
change in geometry in the relaxation method factors in to 
the results. Therefore we suspect that in this case the linear 
methods slightly over-predicting the deflection. 



Table 1 . Maximum displacements of thin and thick cantilever beams 



Geometry 


Mass/spring 


Direct stiffness 


Analytical 


Thin can- 


0.822mm 


0.823mm 


0.823mm 


tilever beam 








Thick can- 


0.538mm 


0.546mm 


N/A 


tilever beam 










(a) 




(b) 

Fig. 4. The deflection of a thick cantilever beam as calculated by 
the direct stiffness method (finite element analysis) (a), and the 
mass/spring method (b). 

To verify the dynamic properties of the simulation, a 
thin cantilever beam of the same dimensions and properties 
as used in the previous section was excited with an impulse 
force at the free end. Damping was turned off, except a trace 
amount of local bond rotational damping (^ = 0.01) to main- 
tain numerical stability. 20,000 data points were collected 
of the z position of the voxel at the free end of the beam, 
corresponding to 0.13 seconds of physical time, or about 60 
oscillations of the lowest-frequency fundamental mode. The 
frequency characteristics of this data were then plotted with 
the analytically calculated natural frequencies of a thin can- 
tilever beam. (Figure [5]) The analytical natural frequencies 
were calculated according to 



G>n — K n 



EI 

fhL 4 



(41) 



where co„ is the natural frequency of mode n in radians 
per second, K n is the standard scaling factor for this mode, 
and m is the mass per unit length of the beam (19). These val- 
ues are overlaid on Figure [5] and the simulated and analytical 



Table 2. The modal frequencies of a simulated thin beam and ana- 
lytically calculated results agree well, with the dynamic mass/spring 
simulation under-predicting slightly due to damping and having a fi- 
nite thickness. 



Mode 


Analvtical 
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2428Hz 


2531Hz 
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6749Hz 
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Fig. 6. The computational speed per voxel drops off slightly as the 
number of voxels in the simulation increases. 
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Fig. 5. The frequency response of a simulated cantilever beam with 
low damping clearly shows modal resonances that agree well with 
analytically calculated values (red overlay lines). 

natural frequencies are tabulated in Table [2] The simulation 
predicts natural frequencies that are slightly lower than those 
predicted by beam theory. This is likely because the 20x1 
aspect ratio of the simulated beam is not quite an ideal thin 
beam, and because there is a small amount of damping in 
the simulation which would tends towards under-predicting 
natural frequencies. 



4 Results 

4.1 Simulation Performance 

Several parameters were explored to characterize the 
performance of the soft body simulator. All results presented 
here assume the simulation is run on a single worker thread 
of CoreI7 CPU at 2.67GHz. As implemented, the simula- 
tion proved very computationally efficient. For a reasonable 
size object of 4000 voxels, 122 complete simluation itera- 
tions were completed per second, or approximately 500,000 
voxel calculations per second. As the number of voxels of 
increases in the object, the total voxels calculated per second 
decreases, but not dramatically. The simulation speed per 
voxel for cubic blocks of various numbers of elements are 
shown in Figure [6] 
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Fig. 7. Damping individual bonds critically (B z = 1.0) lowers the 
noise floor by approximately 7 orders of magnitude compared to the 
undamped case {B z = 0). 

4.2 Effects of local bond damping 

By applying a combination of damping both to the in- 
dividual bonds and to the voxels relative to ground, the so- 
lution converges quickly to steady state with very little nu- 
merical jitter. The addition local damping does not signifi- 
cantly affect the convergence speed, but allows the solution 
to converge to a residual static error approximately 7 orders 
of magnitude lower. By suppressing jitter in this manner, 
both static and dynamic solutions are less susceptible to nu- 
merical instability. 

4.3 Speedup of self collision schemes 

Different combinations of self collision detection meth- 
ods were directly compared using the test geometry shown 
in Figure [8] All materials in this setup were defined with 
a stiffness of IMPa. The red and blue materials each were 
assigned to have a thermal expansion coefficient with mag- 
nitude of 0.02, although one was positive and one negative. 
The temperature of the environment was then sinusoidally 
varied with an amplitude of 30 degrees, which corresponds 
to a 60% expansion and contraction of the red and blue mate- 
rials 180 degrees out of phase. This sets up a periodic colli- 
sion between the extremities that are repeatedly entering and 
exiting the assigned collision horizon. 




Fig. 8. An arbitrary clapper setup to test net iteration rates with dif- 
ferent collision types. The red and blue materials change volume 
sinusoidally 180 degrees out of phase to provide the actuation. 

Table 3. Iteration rates for various collisions detection and handling 
schemes (Higher is better) 



Geometry 


Rate (Iter/sec) 


All+Every 


140.4 


Surf+Every 


140.6 


All+Horizon 


834.4 


Surf+Horizon 


835.4 



Additionally, the incorporating the collision horizon has 
extremely variable speedup up the simulation. In the ex- 
treme case of voxels moving very fast, no speedup is ob- 
served, since the collision horizon may be exceeded every 
timestep. However, this is only possible in extremely fast 
rigid body motion collisions, for which this simulation is not 
intended. On the opposite end of the spectrum, if the object 
is stationary or moving slowly, the collision horizon will take 
many, many timesteps to be exceeded, so almost no collision 
detection calculations will be needed, resulting in dramatic 
acceleration of the simulation. 



5 Demonstrations of simple volumetrically actuated 
mechanisms 

Several demonstration scenes were created to illustrate 
the simulator in action. In the first scene, (Figure [9^-c) an 
actuated beam kicks a ball into a bowling pin. The beam 
has a stiffness 10 times great than the ball, which in turn has 
a stiffness 10 times greater than the bowling pin. The fre- 
quency of the red and blue volumetric actuation was selected 
such that the beam would swing in resonance. The second 
scene (Figure [9]l-f) shows a 2D layer of voxels falling and 
interacting with a fixed sphere as cloth would. In the third 
scene, a quadruped with periodic leg actuation walks forward 
using the nonlinearities of the surface friction with the floor 
as well as the side-to- side resonance of the head swinging 
back and forth. These illustrations demonstrate not just the 
dynamics and large deformation of capabilities of the simu- 
lation, but also the use of volumetric actuation and the robust 
yet efficient collision system. 



The case of comparing the distance from all voxels to 
all other voxel in the structure at every timestep (All+Every) 
was included as a baseline. Comparing only surface vox- 
els to each other at every timestep (Surf+Every) resulted in 
a very minor speedup. This is expected given the geom- 
etry chosen, because the majority of voxels in the struc- 
ture are surface voxels. However, large gains in speed are 
realized when incorporating the collision horizon. Even 
when comparing all voxels to all voxels when recalculation 
is needed (All+Horizon) the simulation as a whole speeds 
up almost 6x. Again, minor acceleration is realized when 
considering only surface voxels with the collision horizon 
(Surf+Horizon). 

It should be stressed that these results are merely repre- 
sentative. In cases with larger numbers of voxels, the bot- 
tleneck of the entire simulation is in collision detection. In 
this case, comparing all voxels to all other voxels is 0(N 2 ) 
(where N is the number of voxels) which dominates the 0(N) 
scaling of calculating forces and updating positions. In this 
case, precompiling a list of surface voxels and only using 
them in collision detection can reduce collision detection 
down to approximately 0(N 1 - 5 ), although this depends on 
the geometry. Of course, if the geometry is thin such that all 
voxels are on the surface, no speedup will be observed. 



6 Conclusion 

We have demonstrated a computationally efficient soft 
body simulator with applications in non-linear material mod- 
eling and dynamic soft object simulation. By virtue of being 
voxel-based, this simulation can accurately model heteroge- 
neous materials with differing stiffnesses and densities in a 
physically accurate environment, even when the materials 
are well-interspersed among each other. This enables mod- 
eling of gradient and composite materials. By incorporating 
a collision horizon that is updated only when needed, self- 
intersection is eliminated with low computational overhead. 
By implementing volumetric actuation, structures can be ac- 
tuated without imposing arbitrary external forces to create 
self-contained mechanisms. 

All code and documentation is freely available at 
www.voxcad.com, including a standalone GUI for editing 
and simulating objects in a real-time interactive environment. 
This simulation opens the door to the design automation of 
a wide variety of non-linear physical structures and mecha- 
nisms that were not possible with previous soft-body physics 
simulation packages. 




(a) (b) 




(d) (e) (f) 




(g) (h) (i) 



Fig. 9. Frames from demonstration scenes show a volumetrically 
actuated flexible beam swung in resonance to kick a soft ball into an 
even softer pin (a-c), A 2D layer of voxels falling under gravity and 
interacting with a fixed sphere (d-f), and a locomoting quadruped (g- 
0- 
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A Simple implementation of the simulator 

A minimal C++ program is provided here to demonstrate the ease of coding a bare-bones dynamic simulation. A 
5mmxl0mmx5mm beam with a material stiffness of lOMPa is cantilevered and a lkN force is applied to the free end. 

CVX_Object Object; //Voxel object 

CVX_Environment Environment; //Environment object 
CVX_Sim Simulator; //Simulator object 

Object . InitializeMatter (0.001 f , 5, 10, 5); //Creates a 5x10x5 workspace of 1mm voxels 
int Matlndex = Object . AddMat( " Material 1 " , 10000000. Of, 0.3f); //Adds a lOMPa material to the palette 
for (int x = 0; x<5; x++){ 
for (int y = 0; y<10; y++){ 
for (int z=0; z<5; z++){ 

Object . SetMat (x , y, z, Matlndex); //Sets each voxel to lOMPa material 

} 

} 

} 

Environment . AddObject(&Obj ect ) ; //Imports the object into the environment 

Environment . AddFixedRegion ( Vec3D (0 , 0, 0), Vec3D(1.0f, O.Olf, l.Of)); //Fixes the -Y plane 
Environment. AddForcedRegion ( Vec3D (0 , 0.99, 0), Vec3D(1.0f, O.Olf, l.Of), Vec3D(0, 0, -1000. Of)); //Adds 
a downward force of lkN to the +Y plane 

Simulator . Import(&Environment ) ; //Imports the environment into the simulator 
Simulator . SetSlowDampZ (0 . 1 3 ) ; //Sets the global damping ratio to an appropriate value 
for(int i=0; i <400; i ++) Simulator . TimeStep () ; //Simulates 400 timesteps 



B Selected function definitions from dynamic voxel simulation classes 

A small subset of functions used in the voxel simulator are documented here. With only these functions a dynamic simulation 
may be created and run. Complete documentation of the source code is available at www.voxcad.com. 

B.l class CVX_Object: 

Describes the geometry and materials of a voxel object. 



void CVX_Object :: InitializeMatter (float iLattice_Dim, int xV, int yV, int zV) 

Initializes voxel object with the specified voxel size and default lattice. A cubic lattice is assumed at the provided inter-voxel 
lattice dimension. 



in 


iLatticeJDim 


The base lattice dimension between adjacent voxels in meters. 


in 


xV 


The number of voxels in the X dimension of the workspace. 


in 


yV 


The number of voxels in the Y dimension of the workspace. 


in 


zV 


The number of voxels in the Z dimension of the workspace. 



int CVX_Object : :AddMat (std: : string Name, double EMod, double PRatio, std: : string *RetMessage = NULL) 

Appends a material to the palette. Returns the index within the palette that this material was created at. This index may 
change if materials are deleted from the palette. All colors and physical properties besides elastic modulus and Poisson's 
ratio are set to defaults. 



in 


Name 


The name of the material to add. If the name is already in use a variant will be generated. 


in 


EMod 


The elatic modulus of the new material in Pascals. 


in 


PRatio 


The poissons ratio of the new material. 


out 


RetMessage 


Pointer to an initialized string. Messages generated in this function will be appended to the string. 



bool CVX_Ob j ect :: SetMat (int x, int y, int z, int Matlndex) [inline] 

Sets a single voxel to the specified material. Returns true if successful. Returns false if provided indices are outside of the 
workspace or the material index is not contained within the current palette. 



in 


X 


Integer X index of voxel location to set. 


in 


y 


Integer Y index of voxel location to set. 


in 


z 


Integer Z index of voxel location to set. 


in 


Matlndex 


Specifies the index within the material palette to set this voxel to. 



B.2 class CVX_Environment : 

Describes the physical environment and boundary conditions for the voxel object. 

void CVX_Environment : : AddOb ject (CVX_Object *p0bjln) [inline] 

Links a voxel object to this environment. Only one voxel object may be linked at a time. 
| in I pObjln | Pointer to an initialized voxel object to link to this simulation. 



void CVXJSnvironment : : AddFixedRegion (Vec3D Location, Vec3D Size) 

Adds a region of voxels to be fixed to ground. All voxels touching this region will be immobile in the simulation. 



in 


Location 


The corner of the region closest to the origin. Specified as a percentage (in X, Y, and Z respectively) of the overall workspace. 
(Location.x, Location.y, and Location.z each have a range of [0.0 to 1.0]). 


in 


Size 


The size of the region. Specified as a percentage (in X, Y, and Z respectively) of the overall workspace. (Size.x, Size.y, and 
Size.z each have a range of [0.0 to 1.0]). 



void CVXJSnvironment : : AddForcedRegion (Vec3D Location, Vec3D Size, Vec3D Force) 

Applies a force to a region of voxels. All voxels touching this region will have an external force applied to them. The 
provided force vector will be divided equally between all voxels. 



in 


Location 


The corner of the region closest to the origin. Specified as a percentage (in X, Y, and Z respectively) of the overall workspace. 
(Location.x, Location.y, and Location.z each have a range of [0.0 to 1.0]). 


in 


Size 


The size of the region. Specified as a percentage (in X, Y, and Z respectively) of the overall workspace. (Size.x, Size.y, and 
Size.z each have a range of [0.0 to 1.0]). 


in 


Force 


The force to be distributed across this region in Newtons. The force is divided equally among all voxels in the forced region. 



B.3 class CVX_Sim: 

Contains the current physical state of the voxel object and advances the simulation. 



void CVX_Sim: : Import (CVX_Environment *pEnvIn = NULL, std: : string *RetMessage = NULL) 

Imports a physical environment into the simulator. The environment should have been previously initialized and linked with 
a single voxel object. This function sets or resets the entire simulation with the new environment. 



in 


pEnvIn 


A pointer to initialized CVX_Environment to import into the simulator. 


out 


RetMessage 


A pointer to initialized string. Output information from the Import function is appended to this string. 



void CVX_Sim: : SetSlowDampZ (double SlowDampIn) [inline] 

Sets the damping ratio that slows downs voxels. When this is non-zero, each voxel is damped (based on its mass and 
stiffness) to ground. Range is [0.0 to 1.0]. Values greater than 1.0 may cause numerical instability. 

| out I SlowDampIn | Damping ratio for damping each voxel relative to ground. | 



void CVX_Sim: : SetCollisionDampZ (double ColDampZIn) [inline] 

Sets the damping ratio for voxels in colliding state. When this is non-zero, each voxel is damped (based on its mass and 
stiffness) according to the relative penetration velocity. Range is [0.0 to 1.0]. Values greater than 1.0 may cause numerical 
instability. 

| out | ColDampZIn | Collision damping ratio for colliding voxels. | 



void CVX_Sim: : SetBondDampZ (double BondDampZIn) [inline] 

Sets the damping ratio for connected voxels. When this is non-zero, each voxel is damped (based on its mass and stiffness) 
according to its relative velocity to the other voxel in each bond. Range is [0.0 to 1.0]. Values greater than 1.0 may cause 
numerical instability. 

| out | BondDampZIn | Damping ratio between each connected voxel. | 



bool CVX_Sim: :TimeStep (std:: string *pRetMessage = NULL) 

Advances the simulation one time step. Given the current state of the simulation (voxel positions and velocities) and infor- 
mation about the current environment, this function advances the simulation by the maximum stable timestep. Returns true 
if the time step was successful, false otherwise. 



out I pRetMessage | Pointer to an initialized string. Messages generated in this function will be appended to the string. 



